function obj = calculate_obj(H, HPs, Ks, beta, gamma, eta, mid_layer_width, layer_connections)
nLayer = length(mid_layer_width);

if nLayer == 1
    iLayer = nLayer;
    gamma_a = gamma{iLayer};
    layer_a = layer_connections{iLayer};
    eta_a = eta{iLayer};
    o1 = 0;
    HPt = HPs{iLayer};
    for i1 = 1 : mid_layer_width(iLayer)
        tmp = HPt{i1} * HPt{i1}';
        tmp = (tmp + tmp')/2;
        gamma_a_b = gamma_a{i1};
        layer_a_b = layer_a{i1};
        for i2 = 1:length(gamma_a_b)
            o1 = o1 + eta_a(i1) * gamma_a_b(i2) * sum(sum( tmp .* Ks(:, :, layer_a_b(i2))));
        end
    end
    HPt = HPs{iLayer};
    o2 = 0;
    tmp = H * H';
    tmp = (tmp + tmp')/2;
    for i1 = 1 : mid_layer_width(iLayer)
        tmp2 = HPt{i1, 1} * HPt{i1, 1}';
        tmp2 = (tmp2 + tmp2')/2;
        o2 = o2 + beta(i1) * sum(sum( tmp .* tmp2));
    end
    obj = o1 + o2;
else
    obj_layers = zeros(nLayer, 1);
    for iLayer = 1:nLayer
        if iLayer == 1
            gamma_a = gamma{iLayer};
            layer_a = layer_connections{iLayer};
            eta_a = eta{iLayer};
            o1 = 0;
            HPt = HPs{iLayer};
            for i1 = 1 : mid_layer_width(iLayer)
                tmp = HPt{i1} * HPt{i1}';
                tmp = (tmp + tmp')/2;
                gamma_a_b = gamma_a{i1};
                layer_a_b = layer_a{i1};
                for i2 = 1:length(gamma_a_b)
                    o1 = o1 + eta_a(i1) * gamma_a_b(i2) * sum(sum( tmp .* Ks(:, :, layer_a_b(i2))));
                end
            end
            obj_layers(iLayer) = o1;
            
        elseif iLayer == nLayer
            HPt = HPs{iLayer};
            o1 = 0;
            tmp = H * H';
            tmp = (tmp + tmp')/2;
            for i1 = 1 : mid_layer_width(iLayer)
                tmp2 = HPt{i1, 1} * HPt{i1, 1}';
                tmp2 = (tmp2 + tmp2')/2;
                o1 = o1 + beta(i1) * sum(sum( tmp .* tmp2));
            end
            obj_layers(iLayer) = o1;
            
        else
            gamma_a = gamma{iLayer};
            layer_a = layer_connections{iLayer};
            eta_a = eta{iLayer};
            o1 = 0;
            HPt0 = HPs{iLayer - 1};
            HPt = HPs{iLayer};
            for i1 = 1 : mid_layer_width(iLayer)
                tmp = HPt{i1} * HPt{i1}';
                tmp = (tmp + tmp')/2;
                gamma_a_b = gamma_a{i1};
                layer_a_b = layer_a{i1};
                for i2 = 1:length(gamma_a_b)
                    tmp2 = HPt0{layer_a_b(i2), 1} * HPt0{layer_a_b(i2), 1}';
                    tmp2 = (tmp2 + tmp2')/2;
                    o1 = o1 + eta_a(i1) * gamma_a_b(i2) * sum(sum( tmp .* tmp2 ));
                end
            end
            obj_layers(iLayer) = o1;
        end
    end
    obj = sum(obj_layers);
end
end